##################################################
#
# PLOTS OF SOCIAL VICE (ALCOHOL AND TOBACCO) MODELS FOR
# OIL, ISLAM, and WOMEN PROJECT
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Created 27 February 2016
#
##################################################

##################################################
# R Housekeeping
##################################################

# Remove and empty workspace
rm(list=ls())

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Packages
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

library(foreign)
library(sandwich)
library(apsrtable)
library(tonymisc)  # includes summaryR()
library(MASS) # for polr() which runs our ologit

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Functions
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 


# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Code",sep=""))
source("statifyFN.R")
source("multipleOLSFN.R")
source("genModelsFN.R")
source("plotlineFN.R")
source("plotprepFN.R")
source("modelDefinitions.R")
source("tableOLSFN.R")

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load DataISQFINAL
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd(paste(mywd,"/Data",sep=""))

data <- read.dta("NewDVModelsXSforISQFINAL2015.dta")

#################################################
#
# Models of Oil Income and Reproductive Issues
# 
#################################################


out.malesmoking <- multipleOLS(formulas=genModels(c("malesmoking"),
                                            c(p142.09)),
				          	modelnames=c(names(p142.09)))

out.femalesmoking <- multipleOLS(formulas=genModels(c("femsmoking"),
                                             c(p142.09)),
                          modelnames=c(names(p142.09)))

out.alcohol <- multipleOLS(formulas=genModels(c("totalalcohollitres"),
                                             c(p140.04,p142.04)),
                          modelnames=c(LETTERS[1:7]))



## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
# Plotting Models
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
vars.to.plot <- c("oil_gas_valuePOPred","OIPC0009_1k",
                  "oil_gas_pc_09","oil_gas_pc_04")
				
plot.out.malesmoking <- do.call(rbind.data.frame,
                            lapply(out.malesmoking,plotprep,vars.to.plot))
plot.out.femalesmoking <- do.call(rbind.data.frame,
                           lapply(out.femalesmoking,plotprep,vars.to.plot))
plot.out.alcohol <- do.call(rbind.data.frame,
                            lapply(out.alcohol,plotprep,vars.to.plot))

plot.out.malesmoking$Index <- 1:nrow(plot.out.malesmoking)
plot.out.femalesmoking$Index <- 1:nrow(plot.out.femalesmoking)
plot.out.alcohol$Index <- 1:nrow(plot.out.alcohol)

plot.out.malesmoking$bg <- ifelse(abs(plot.out.malesmoking$Estimate/plot.out.malesmoking$SE)<=1.65, "white",
                                  ifelse(abs(plot.out.malesmoking$Estimate/plot.out.malesmoking$SE)>=1.96,"black",
                                         "gray"))
plot.out.femalesmoking$bg <- ifelse(abs(plot.out.femalesmoking$Estimate/plot.out.femalesmoking$SE)<=1.65, "white",
                                  ifelse(abs(plot.out.femalesmoking$Estimate/plot.out.femalesmoking$SE)>=1.96,"black",
                                         "gray"))
plot.out.alcohol$bg <- ifelse(abs(plot.out.alcohol$Estimate/plot.out.alcohol$SE)<=1.65, "white",
                                    ifelse(abs(plot.out.alcohol$Estimate/plot.out.alcohol$SE)>=1.96,"black",
                                           "gray"))

# pch = 21
# col = black

setwd("../Drafts/Charts")

#pdf(file="ISQFINALMainBodyFigure6.pdf",height=5,width=7)
tiff(file="ISQFINALMainBodyFigure6.tiff",
			height=5,
			width=7,
			res=300,
			units="in")
			
par(mfrow=c(3,1),oma=c(0,0,3,1),mar=c(3,3,1,1))

plot(3,
     pch="",
     ylim=c(1.15*min(plot.out.alcohol$Y0),
            abs(1.15*max(plot.out.alcohol$Y1))), 
     xlim=c(.5,max(plot.out.alcohol$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Alcohol Consumption",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.alcohol$Y0,
         y1=plot.out.alcohol$Y1,
         x0=plot.out.alcohol$Index,
         x1=plot.out.alcohol$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.alcohol$Index,
       y=plot.out.alcohol$Estimate,
       col="black",
       pch=21,
       bg=plot.out.alcohol$bg,
       cex=1.5)
axis(side=1,at=plot.out.alcohol$Index,
     labels=plot.out.alcohol[1:nrow(plot.out.alcohol),]$Model,
     cex.axis=1.25, tick=TRUE)

plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.malesmoking$Y0),
            abs(1.15*max(plot.out.malesmoking$Y1))), 
     xlim=c(.5,max(plot.out.malesmoking$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Male Smoking",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.malesmoking$Y0,
		y1=plot.out.malesmoking$Y1,
		x0=plot.out.malesmoking$Index,
		x1=plot.out.malesmoking$Index,
		col="black",
		lwd=1,
		lty=1)
points(x=plot.out.malesmoking$Index,
		y=plot.out.malesmoking$Estimate,
		col="black",
		pch=21,
		bg=plot.out.malesmoking$bg,
		cex=1.5)
axis(side=1,at=plot.out.malesmoking$Index,
		labels=plot.out.malesmoking[1:nrow(plot.out.malesmoking),]$Model,
		cex.axis=1.25, tick=TRUE)
		

### And Female Smoking
plot(3,
     pch="",
     ylim=c(1.15*min(plot.out.femalesmoking$Y0),
            abs(1.15*max(plot.out.femalesmoking$Y1))), 
     xlim=c(.5,max(plot.out.femalesmoking$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Estimated Relationship With Female Smoking",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.femalesmoking$Y0,
         y1=plot.out.femalesmoking$Y1,
         x0=plot.out.femalesmoking$Index,
         x1=plot.out.femalesmoking$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.femalesmoking$Index,
       y=plot.out.femalesmoking$Estimate,
       col="black",
       pch=21,
       bg=plot.out.femalesmoking$bg,
       cex=1.5)
axis(side=1,at=plot.out.femalesmoking$Index,
     labels=plot.out.femalesmoking[1:nrow(plot.out.femalesmoking),]$Model,
     cex.axis=1.25, tick=TRUE)


# Title for whole thing
mtext("Estimating Relationship Between Oil Income And Other Social Behaviors", 
      side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()



###